{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matrix Algebra,\n",
    "\n",
    "## A first introduction from the computational standpoint\n",
    "\n",
    "This notebook is part of the [Virgilio project](https://github.com/clone95/Virgilio) and serves as a first introduction to matrices and how they work. You are welcome to read this at any level of knowledge, but I would recommend it to people who\n",
    "\n",
    " - want to see how matrices work\n",
    " - want to implement a basic Python library to deal with matrices\n",
    " - have already dealt with matrices but need a refresher\n",
    " \n",
    "Conversely, I would not classify this as the best resource for you if\n",
    "\n",
    " - you want to learn the deep insights of how matrices work and why they work the way they do\n",
    " - you want to get a nice knowledge of linear algebra: matrices are just a part of linear algebra, not the whole thing!\n",
    " \n",
    "This guide is broken down in several sections, as you can see from the following index:\n",
    "\n",
    " - [Further reading](#Further-reading)\n",
    " - [What is a matrix?](#What-is-a-matrix?)\n",
    "   - [Rows and columns](#Rows-and-columns)\n",
    " - [Arithmetics with matrices](#Arithmetics-with-matrices)\n",
    "   - [Addition and subtraction](#Addition-and-subtraction)\n",
    "   - [Multiplication](#Multiplication)\n",
    "     - [Vector dot product](#Vector-dot-product)\n",
    "     - [Back to multiplication](#Back-to-multiplication)\n",
    "     - [Identity matrix](#Identity-matrix)\n",
    "     - [Scalar multiplication](#Scalar-multiplication)\n",
    "   - [Division..???](#Division..???)\n",
    " - [Inverse](#Inverse)\n",
    "   - [Determinant](#Determinant)\n",
    " - [Transpose](#Transpose)\n",
    " - [Linear systems](#Linear-systems)\n",
    "   - [Gaussian elimination](#Gaussian-elimination)\n",
    " - [Inverses revisited](#Inverses-revisited)\n",
    " - [Trace](#Trace)\n",
    " - [Final Matrix implementation](#Final-Matrix-implementation)\n",
    " - [Contacts](#Contacts)\n",
    " \n",
    "If at any point you find a problem with the code, please refer to the [Contacts](#Contacts) section and write to the author or directly address the problem by [submitting a pull request](https://github.com/clone95/Virgilio/pulls) or [opening an issue](https://github.com/clone95/Virgilio/issues) on the [GitHub page of project Virgilio](https://github.com/clone95/Virgilio)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A forword on the length of this notebook\n",
    "\n",
    "If you look at the scroll bar, this is one heck of a lengthy notebook! Don't worry, really: there is basically a _lot_ of duplicated code from the `Matrix` class we will be creating, because whenever we write a new function for it I just include the whole class again (this way, any particular cell of code can always be ran independently and you don't have to evaluate all the previous cells for it to work) and because of this, the notebook grew a bit bigger than needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further reading\n",
    "\n",
    "Usually the _Further reading_ section comes in the end... I decided to include it in the beginning because you may find it better to first read something more general or beginner friendly and then read this, as the focus of this document will be implementing the usual matrix operations.\n",
    "\n",
    "If you just want to know what a matrix is, you can check [Wikipedia's article](https://en.wikipedia.org/wiki/Matrix_(mathematics)) on matrices. If you are interested in the definition of matrix from a linear algebra standpoint, you can [check it](http://mathworld.wolfram.com/Matrix.html) on WolframAlpha's website; if you scroll down there is a **SEE ALSO** section that will point you to lots of interesting things.\n",
    "\n",
    "If you are confused about linear algebra and working with matrices, at the time of writing (28/04/2019) the top 3 answers of [this MO question](https://mathoverflow.net/questions/11669/what-is-the-difference-between-matrix-theory-and-linear-algebra) are quite enlightening. In particular, [this answer](https://mathoverflow.net/a/11679) contains the following excerpt:\n",
    " > [...] A matrix is just a list of numbers, and you're allowed to add and multiply matrices by combining those numbers in a certain way. [...]  You can do lots of interesting numerical things with matrices, and these interesting numerical things are very important because matrices show up a lot in engineering and the sciences.\n",
    " > In linear algebra, however, you instead talk about linear transformations, which are not (I cannot emphasize this enough) a list of numbers, although sometimes it is convenient to use a particular matrix to write down a linear transformation.\n",
    " \n",
    "[The third answer](https://mathoverflow.net/a/19884) contains a quotation from a mathematician, regarding the fact that matrices are usually used as a tool for linear algebra:\n",
    " > \"There is hardly any theory which is more elementary [than linear algebra], in spite of the fact that generations of professors and textbook writers have obscured its simplicity by preposterous calculations with matrices.\"\n",
    " \n",
    "And [this answer](https://mathoverflow.net/a/19923) includes yet another quotation, where a mathematician confesses that even though matrices are just tools of linear algebra, they are quite useful:\n",
    " > \"We share a philosophy about linear algebra: we think basis-free, we write basis-free, but when the chips are down we close the office door and compute with matrices like fury.\"\n",
    " \n",
    "Either way, even if you do or do not understand what is being exposed on the previous excerpts, [KhanAcademy](https://www.khanacademy.org/math/linear-algebra) has a short course on Linear Algebra, which is coupled with an intro to matrices. (For those of you who don't know KhanAcademy, I highly recommend it; it is one of my favourite e-learning resources.) The Youtube channel [3blue1brown](https://www.youtube.com/channel/UCYO_jab_esuFRV4b17AJtAw) also has a [nice series on Linear Algebra](https://www.3blue1brown.com/essence-of-linear-algebra-page) with every single concept explained visually.\n",
    "\n",
    "Lastly, Reddit's [/r/learnmath](https://www.reddit.com/r/learnmath/) subreddit has a post with a [very nice collection](https://www.reddit.com/r/learnmath/comments/8p922p/list_of_websites_ebooks_downloads_etc_for_mobile/) of resources, among them resources for Linear Algebra."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is a matrix?\n",
    "\n",
    "Bluntly put, a matrix is a table of numbers, arranged in rows and columns:\n",
    "\n",
    "$$\\begin{bmatrix} 1 & 2 \\\\ 3 & 4 \\end{bmatrix}$$\n",
    "\n",
    "You can also see them represented as\n",
    "\n",
    "$$\\begin{pmatrix} 1 & 2 \\\\ 3 & 4 \\end{pmatrix}$$\n",
    "\n",
    "and using parenthesis or brackets is up to your personal taste, but please refrain from writing matrices as\n",
    "\n",
    "$$\\begin{vmatrix} 1 & 2 \\\\ 3 & 4 \\end{vmatrix}$$\n",
    "\n",
    "as those vertical bars have a special meaning.\n",
    "\n",
    "The number of _rows_ of a matrix is the number of \"horizontal lines\" the matrix has and the number of _columns_ of a matrix is the number of \"vertical lines\" the matrix has. If a matrix has $r$ rows and $c$ columns we usually write that the matrix is a $r \\times c$ matrix, which we read as _$r$ by $c$_ matrix. In mathematical notation, it is standard to represent matrices as uppercase letters like $A$, $B$, $M$.\n",
    "\n",
    "We will start our implementation now, and we will define a matrix as a Python object that takes its number of rows and columns as input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Matrix (3x5)>\n",
      "3\n",
      "[<M (3x5)>, <M (3x5)>, <M (1x10)>]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "matrix = Matrix(3, 5)\n",
    "print(matrix)\n",
    "print(matrix.nrows())\n",
    "list_of_matrices = [matrix, matrix, Matrix(1, 10)]\n",
    "print(list_of_matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we also saw that matrices are supposed to hold numbers, so we will add that functionality as well. We will use _a list of lists_ to save the values internally, and we will say that each list inside the `values` variable is a row of the matrix, so that\n",
    "\n",
    "`values = [\n",
    "    [1, 2],\n",
    "    [3, 4]\n",
    "]`\n",
    "\n",
    "matches the matrix\n",
    "\n",
    "$$\\begin{bmatrix} 1 & 2 \\\\ 3 & 4 \\end{bmatrix}$$\n",
    "\n",
    "So our implementation will take an extra optional argument `values` which gives the initial matrix values to consider. If nothing is provided, we will initialize everything with $0$.\n",
    "\n",
    "So this new implementation contains:\n",
    " 1. an improved `__init__` function on lines 2-25\n",
    "\n",
    "(to see the line numbering, select a cell and hit \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "values = [[1,2], [3,4]]\n",
    "matrix = Matrix(2, 2, values)\n",
    "# this throws an error because the sizes are wrong\n",
    "#matrix = Matrix(10, 34, values)\n",
    "# giving no values still works:\n",
    "matrix = Matrix(2, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rows and columns\n",
    "\n",
    "Like I mentioned previously, a matrix $A$ always has a certain number of rows and a certain number of columns. As we will start to see shortly, it is often convenient to refer to specific rows/columns of our matrix. For the rows, we start to count them from the top and for the columns, we start to count them from the left. So if your matrix is\n",
    "\n",
    "$$\\begin{bmatrix} 1 & 2 & 3 \\\\ 4 & 5 & 6 \\\\ 7 & 8 & 9\\end{bmatrix}$$\n",
    "\n",
    "then its second row is\n",
    "\n",
    "$$\\begin{bmatrix} 4 & 5 & 6 \\end{bmatrix}$$\n",
    "\n",
    "and its third column is\n",
    "\n",
    "$$\\begin{bmatrix} 3 \\\\ 6 \\\\ 9 \\end{bmatrix}$$\n",
    "\n",
    "Notice that, in theory, the orientation matters! That is,\n",
    "\n",
    "$$\\begin{bmatrix} 3 \\\\ 6 \\\\ 9 \\end{bmatrix} \\neq \\begin{bmatrix} 3 & 6 & 9 \\end{bmatrix}$$\n",
    "\n",
    "Strictly speaking, a matrix with only one row is called a row vector and a matrix with only one column is called a column vector. Both of them are vectors. There are times when distinguishing between row/column vectors makes no difference and there are times when distinguishing the two makes **all** the difference!\n",
    "\n",
    "We will now expand our implementation to include a way of retrieving a specific row/column and we will also include functions to get a list of all rows and a list of all columns. **We will follow the mathematical convention that the first row is row 1** and so the indexing of rows, in the mathematical sense, and of lists in the implementation, has to be carefully handled.\n",
    "\n",
    "This new implementation contains:\n",
    " - 27:32, a helper function `from_values` that allows creating a Matrix by just providing the initial list of values\n",
    " - 38:50, a function `get_row` to return a specific row of the matrix\n",
    " - 52:54, a function `get_rows` to return all rows of the matrix\n",
    " - 60:72 and 74:76, the corresponding `get_col` and `get_cols` functions for columns\n",
    "\n",
    "(to see the line numbering, select a cell and hit \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7, 8, 9]\n",
      "[1, 4, 7]\n",
      "[[1, 4, 7], [2, 5, 8], [3, 6, 9]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "values = [[1,2,3],[4,5,6],[7,8,9]]\n",
    "matrix = Matrix.from_values(values)\n",
    "print(matrix.get_row(3))\n",
    "print(matrix.get_col(1))\n",
    "# this raises an error\n",
    "#print(matrix.get_row(0))\n",
    "# and so does this\n",
    "#print(matrix.get_row(10))\n",
    "print(matrix.get_cols())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arithmetics with matrices\n",
    "\n",
    "Now that we have defined our matrices, it is time to do operations with them! The most basic operation we can do is add two matrices together.\n",
    "\n",
    "### Addition and subtraction\n",
    "\n",
    "And how would that work? When we are adding two matrices together, we add the values on the corresponding positions, so that\n",
    "\n",
    "$$\\begin{bmatrix} 1 & 2 \\\\ 3 & 4\\end{bmatrix} + \\begin{bmatrix} 1 & -1 \\\\ 3 & 73\\end{bmatrix} = \\begin{bmatrix} 1+1 & 2 + (-1) \\\\ 3 + 3 & 4 + 73\\end{bmatrix} = \\begin{bmatrix} 2 & 1 \\\\ 6 & 77\\end{bmatrix}$$\n",
    "\n",
    "because of this, it makes no sense to try to add two matrices that do not have the same _shape_. When we talk about shape, we are talking the number of rows and number of columns. So if your matrix is $n \\times m$ (remember, this means the matrix has $n$ rows and $m$ columns) you can only add it to another $n \\times m$ matrix. As expected, the same goes with subtraction.\n",
    "\n",
    "Now that we know about this, we want to be able to extend our implementation of the matrices. In Python, when we want to implement addition and subtraction, we can implement the [magic methods](https://rszalski.github.io/magicmethods/) `__add__` and `__sub__`. These methods will return a *new* matrix with the right values.\n",
    "\n",
    "We will now implement these two operations. Our new implementation contains:\n",
    " 1. A function `shape` on lines 34-36 that return the shape of a matrix, to check for compatibility\n",
    " 2. The implementation for the function `__add__` on lines 82-96, that adds two Matrix objects together\n",
    " 3. The implementation for the function `__sub__` on lines 98-110, that subtracts one Matrix from another Matrix\n",
    "\n",
    "(to see the line numbering, select a code cell and press \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 2)\n",
      "(10, 20)\n",
      "[[2, 1], [6, 77]]\n",
      "[[0, 3], [0, -69]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "m1 = Matrix.from_values([[1,2],[3,4]])\n",
    "m2 = Matrix.from_values([[1,-1],[3,73]])\n",
    "oddly_shaped = Matrix(10, 20)\n",
    "\n",
    "print(m1.shape())\n",
    "print(oddly_shaped.shape())\n",
    "m = m1 + m2\n",
    "# we shouldn't access _values directly; we'll fix this soon enough\n",
    "print(m._values)\n",
    "m = m1 - m2\n",
    "print(m._values)\n",
    "# this throws an error:\n",
    "#print(m1 + oddly_shaped)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiplication\n",
    "\n",
    "One can also multiply matrices, but that is trickier than it might seem. Can you try to guess what is\n",
    "\n",
    "$$\\begin{bmatrix} 1 & 2 \\\\ 2 & 1 \\end{bmatrix} \\times \\begin{bmatrix} 3 & 3 \\\\ 3 & 3 \\end{bmatrix} =\\ ?$$\n",
    "\n",
    "Well, if you guessed\n",
    "\n",
    "$$\\begin{bmatrix}3 & 6 \\\\ 6 & 3\\end{bmatrix}$$\n",
    "\n",
    "then your guess makes a lot of sense but is wrong; in fact, the right result is\n",
    "\n",
    "$$\\begin{bmatrix} 9 & 9 \\\\ 9 & 9\\end{bmatrix}$$\n",
    "\n",
    "Which looks quite funny, if you've never seen matrix multiplication! Doesn't it? It is. Now I just have to tell you how we get there; but first, let us revisit what the vector dot product is, because that will make it much easier for me to explain matrix multiplication.\n",
    "\n",
    "#### Vector dot product\n",
    "\n",
    "Let us say that $v$ is a vector, something like $v = (1, 3, -2)$; also, let me take $u$ as another vector, say $u = (\\frac12, 1, -1)$. For those of you who are not familiar with it, the dot product of $u$ and $v$, usually written as $u \\cdot v = (\\frac12, 1, -1) \\cdot (1, 3, -2)$ is the number $\\frac12 \\times 1 + 1 \\times 3 + (-1)\\times(-2) = 5.5$\n",
    "\n",
    "The way it works is, we multiply the corresponding values together and then add everything. In mathematical notation, if we have two vectors $v = (v_1, v_2, \\cdots, v_n)$ and $u = (u_1, u_2, \\cdots, u_n)$, then the dot product $u\\cdot v$ is\n",
    "\n",
    "$$u \\cdot v = \\sum_{i = 1}^n u_iv_i$$\n",
    "\n",
    "From the definition it should be clear that $u \\cdot v = v\\cdot u$.\n",
    "\n",
    "Let us implement the dot product as a separate function really quickly; the function `dot_product(u, v)` should accept two vectors and compute their dot product (for us, vectors will be lists):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.5\n",
      "5.5\n"
     ]
    }
   ],
   "source": [
    "def dot_product(u, v):\n",
    "    acc = 0\n",
    "    for ui, vi in zip(u, v):\n",
    "        acc += ui*vi\n",
    "    return acc\n",
    "\n",
    "u = [0.5, 1, -1]\n",
    "v = [1, 3, -2]\n",
    "print(dot_product(u, v))\n",
    "print(dot_product(v, u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the definition, we also get that the dot product only makes sense between two vectors of the same size, i.e. that have the same dimensions! So if $u$ is a vector with $n$ coordinates and if $v$ is a vector with $m$ coordinates, we can only compute $u \\cdot v$ if $n = m$.\n",
    "\n",
    "#### Back to multiplication\n",
    "\n",
    "How can we relate this with multiplication? Well, let us just agree that if I have two vectors $u$ and $v$, then their dot product $u\\cdot v$ is the same as the *matrix* multiplication $u \\times v$ **IF** we write $u$ as a row vector and $v$ as a column vector.\n",
    "\n",
    "That is,\n",
    "\n",
    "$$5.5 = \\frac12 \\times 1 + 1\\times 3 + (-1)\\times(-2) = \\begin{bmatrix}\\frac12 & 1 & -1\\end{bmatrix} \\times \\begin{bmatrix} 1 \\\\ 3 \\\\-2 \\end{bmatrix}$$\n",
    "\n",
    "So, in order to understand what the matrix multiplication $A \\times B$ is, we just see $A$ as a bunch of row vectors stacked on top of eachother and we see $B$ as a bunch of column vectors stacked side-by-side. What I mean is, instead of seeing the matrix multiplication as this:\n",
    "\n",
    "![Two matrices being multiplied](imgbin/matrix_multiplication.png)\n",
    "\n",
    "see the multiplication as the several vectors being combined:\n",
    "\n",
    "![The matrices were split into vectors; the left matrix was split into row vectors and the right matrix was split into column vectors](imgbin/matrix_multiplication_broken.png)\n",
    "\n",
    "So each element in the resulting matrix is the dot product between a row of the left matrix and a column of the right matrix. For example, how can we find this specific value:\n",
    "\n",
    "$$A \\times B = \\begin{bmatrix} 1 & 2 & 3 \\\\ 4 & 5 & 6 \\\\ 7 & 8 & 9\\end{bmatrix} \\times \\begin{bmatrix} 10 & 11 & 12 \\\\ 13 & 14 & 15 \\\\ 16 & 17 & 18 \\end{bmatrix} = \\begin{bmatrix} * & * & * \\\\ * & * & ? \\\\ * & * & * \\end{bmatrix}$$\n",
    "\n",
    "That value is the dot product between a row of $A$ and a column of $B$. Because $?$ is in the second row, we are going to take the second row from $A$, and because $?$ is in the third column, we will take the third column of $B$. Then we do the dot product of those and get\n",
    "\n",
    "$$? = 4\\times12 + 5\\times15 + 6\\times18 = 48 + 75 + 108 = 231 $$\n",
    "\n",
    "so that\n",
    "\n",
    "$$A \\times B = \\begin{bmatrix} 1 & 2 & 3 \\\\ 4 & 5 & 6 \\\\ 7 & 8 & 9\\end{bmatrix} \\times \\begin{bmatrix} 10 & 11 & 12 \\\\ 13 & 14 & 15 \\\\ 16 & 17 & 18 \\end{bmatrix} = \\begin{bmatrix} * & * & * \\\\ * & * & 231 \\\\ * & * & * \\end{bmatrix}$$\n",
    "\n",
    "One final note: the matrices do not need to have the same shape! The only thing that must happen is that we need to take the dot product of rows from the left with columns from the right, that is, rows on the left should have the same length of columns on the right; that means, the number of columns of the left matrix should match the number of rows of the right matrix! If that is the case, then we can take the matrix multiplication. If $A$ is a $n \\times m_A$ matrix and if $B$ is an $m_B \\times p$ matrix, then $A\\times B$ only makes sense when $m_A = m_B$, and in that case the result is a matrix with shape $n \\times p$.\n",
    "\n",
    "With this characterization of matrix multiplication, we can already implement it; below, our `Matrix` implementation added:\n",
    "  1. On lines 112-118, the `dot_product` function\n",
    "  2. On lines 120-134, the `__mul__` magic function, the magic function for multiplication\n",
    "  \n",
    "(to see the line numbering select a cell and then press \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[84, 90, 96]\n",
      "[201, 216, 231]\n",
      "[318, 342, 366]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies two matrices\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Expected a matrix\")\n",
    "        elif self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        # this gives the shape of the result:\n",
    "        nrows = self.nrows()\n",
    "        ncols = other.ncols()\n",
    "        values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "        for r, row in enumerate(self.get_rows()):\n",
    "            for c, col in enumerate(other.get_cols()):\n",
    "                values[r][c] = Matrix.dot_product(row, col)\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "A = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
    "B = [[10, 11, 12], [13, 14, 15], [16, 17, 18]]\n",
    "A = Matrix.from_values(A)\n",
    "B = Matrix.from_values(B)\n",
    "result = A*B\n",
    "for row in result.get_rows():\n",
    "    print(row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that, in general, people will define matrix multiplication directly with its formula, without motivating it with the dot product. If that is the case, you might find something like this written:\n",
    "\n",
    " > If $A$ is an $n \\times m$ matrix and $B$ is an $m \\times p$ matrix, then the product $A \\times B$ is the $n \\times p$ matrix with $$(A\\times B)_{i,j} = \\sum_{k=1}^m a_{i,k}b_{k,j}$$\n",
    " \n",
    "where the $M_{i,j}$ notation means the value of matrix $M$ that is in the $i$-th row, $j$-th column.\n",
    "\n",
    "#### Identity matrix\n",
    "\n",
    "Remember that $1\\times x = x \\times 1$. Is there a similar thing for matrices? Is there some special matrix $I$ that behaves, with matrices, like $1$ behaves with numbers? Is there a matrix $I$ such that\n",
    "\n",
    "$$I \\times M = M \\times I = M,\\quad \\forall M$$\n",
    "\n",
    "Well, we just saw that in order to be able to multiply matrices, their shapes need to match in some way. If we want to be able to do $I \\times M$ and $M \\times I$ then we must have $I$ and $M$ square matrices of the same size; otherwise one of the multiplications won't work. But if we fix the size of $M$ to be $n\\times n$, there is a special matrix $I_n$ such that\n",
    "\n",
    "$$I_n \\times M = M \\times I_n = M,\\quad \\forall M$$\n",
    "\n",
    "It is the matrix that has $0$s everywhere, except in the positions where the index of the row and the column are the same (we call these positions the diagonal of the matrix). Here are $I_n$ for $n \\in \\{1,2,3\\}$:\n",
    "\n",
    "$$I_1 = \\begin{bmatrix} 1 \\end{bmatrix}$$\n",
    "\n",
    "$$I_2 = \\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix}$$\n",
    "\n",
    "$$I_3 = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix}$$\n",
    "\n",
    "We can include a helper method in our `Matrix` object to be able to generate the identity matrix $I_n$ for a given $n$. And while we are at it, recall we saw that the identity matrix is a matrix with a very special format: $0$s everywhere except in the diagonal. (For matrices where this is the case, we call them diagonal matrices and we can also represent them with $D = diag(d_1, d_2, \\cdots, d_n)$ where $d_i = D_{i,i}$.) We will instead create a function that can generate the matrix $D$ from its diagonal values $(d_1, \\cdots, d_n)$ and use that function to create our identity matrix with $(1, \\cdots, 1)$.\n",
    "\n",
    "Our new implementation now contains,\n",
    " 1. on lines 34-40, a function `from_diagonal` that generates a diagonal matrix with the given values in the diagonal\n",
    " 2. on lines 42-45, a function `eye` that returns the identity matrix I (the letter I is read as \"eye\") of the given dimension\n",
    " 3. on lines 157-159, a function `full_str` that returns a string with a printable version of the values of the matrix\n",
    " \n",
    "(to see the line numbering select a cell and then press \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1, 0, 0]\n",
      " [0, 1, 0]\n",
      " [0, 0, 1]]\n",
      "[[1, 0, 0]\n",
      " [0, 4, 0]\n",
      " [0, 0, 7]]\n",
      "[[1, 0, 0]\n",
      " [0, 4, 0]\n",
      " [0, 0, 7]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies two matrices\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Expected a matrix\")\n",
    "        elif self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        # this gives the shape of the result:\n",
    "        nrows = self.nrows()\n",
    "        ncols = other.ncols()\n",
    "        values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "        for r, row in enumerate(self.get_rows()):\n",
    "            for c, col in enumerate(other.get_cols()):\n",
    "                values[r][c] = Matrix.dot_product(row, col)\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\"\n",
    "    \n",
    "eye = Matrix.eye(3)\n",
    "print(eye.full_str())\n",
    "other = Matrix.from_diagonal([1, 4, 7])\n",
    "print((eye*other).full_str())\n",
    "print((other*eye).full_str())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Scalar multiplication\n",
    "\n",
    "Besides multiplying matrices with other matrices, you can also multiply a matrix with a scalar value. That is, if we have a real or complex number, we can multiply any matrix we want with it. If $a$ is a real/complex scalar and $M$ is a matrix, where $(M)_{i,j} = m_{i,j}$ then\n",
    "\n",
    "$$(aM)_{i,j} = am_{i,j}$$\n",
    "\n",
    "The expression above reads: the value in the $i$-th row, $j$-th column of the product $aM$ is $a$ times whatever was in the $i$-th row, $j$-th column of $M$.\n",
    "\n",
    "To implement this operation we just need to change the `__mul__` implementation we have and add the `__rmul__` magic method. For this, we\n",
    " 1. changed the `__mul__` implementation on lines 133-155 to handle scalars\n",
    " 2. added the magic method `__rmul__` on lines 155-158 to handle the cases where the matrix is on the right of a `*`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[60, 20]\n",
      " [30, 615]]\n",
      "[[0.0, 0.0]\n",
      " [0.0, 0.0]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies a matrix with another matrix with a scalar\"\"\"\n",
    "        if not isinstance(other, (Matrix, int, float, complex)):\n",
    "            raise TypeError(\"Cannot multiply Matrix with {}\".format(type(other)))\n",
    "        elif isinstance(other, Matrix) and self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        \n",
    "        if isinstance(other, Matrix):\n",
    "            # this gives the shape of the result:\n",
    "            nrows = self.nrows()\n",
    "            ncols = other.ncols()\n",
    "            values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "            for r, row in enumerate(self.get_rows()):\n",
    "                for c, col in enumerate(other.get_cols()):\n",
    "                    values[r][c] = Matrix.dot_product(row, col)\n",
    "            return Matrix.from_values(values)\n",
    "        else:\n",
    "            return Matrix.from_values( [\n",
    "                [value*other for value in row] for row in self.get_rows()\n",
    "            ])\n",
    "        \n",
    "    def __rmul__(self, other):\n",
    "        \"\"\"This is called when self was on the right-hand side of the * operator\n",
    "            as in 5*Matrix or 4.2*Matrix\"\"\"\n",
    "        return self.__mul__(other)\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\"\n",
    "    \n",
    "m = Matrix.from_values([[12,4],[6, 123]])\n",
    "print((m*5).full_str())\n",
    "print((0.0*m).full_str())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Division..???\n",
    "\n",
    "Now that we know how to add, subtract and multiply matrices, we are left with learning how to divide matrices... Except that such a thing **doesn't exist**! How do you want to define matrix division? We are used to having division _undo_ whatever multitplication does, right? If $y \\neq 0$, then we have\n",
    "\n",
    "$$(x \\times y) \\div y = x$$\n",
    "\n",
    "which we usually write\n",
    "\n",
    "$$\\frac{x \\times y}{y} = x$$\n",
    "\n",
    "So we could hope that matrix _\"division\"_ also works that way, but we don't really know how to _\"divide\"_ a matrix by another, i.e. this is not defined: $A \\div B$ for general matrices. This only makes sense if $B$ is a $1 \\times 1$ matrix, in which case $A \\div B = A \\times \\frac{1}{B}$. Instead, we talk about [matrix inverses](#Inverse)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inverse\n",
    "\n",
    "If $A$ is a matrix, then we say the **inverse of $A$** is the matrix $M$ such that $A \\times M = M \\times A = I$, where $I$ is the identity matrix.\n",
    "\n",
    "From this _definition_, we already get that only square matrices can have an inverse, because the inverse of a matrix $A$, which we usually write $A^{-1}$, should be multipliable by $A$ on both sides. Notice that the notation $A^{-1}$ is similar to the notation we use for real numbers. The inverse of a number $x \\neq 0$ is $\\frac{1}{x}$ which can also be written $x^{-1}$. Also, $x \\times x^{-1} = x^{-1} \\times x = 1$, so that kind of resembles what we wrote for matrices: $A \\times A^{-1} = A^{-1} \\times A = I$.\n",
    "\n",
    "Whenever we talk about matrix inverses, we also have to talk about the [determinant](#Determinant) of a matrix. Notice how I've written that $x \\neq 0$ has an inverse $x^{-1}$? Well, for matrices it is not enough to say that the matrix is not the matrix filled with zeroes. So, in fact, this matrix does not have an inverse:\n",
    "\n",
    "$$\\begin{bmatrix} 0 & 0 \\\\ 0 & 0 \\end{bmatrix}$$\n",
    "\n",
    "But this one does not have an inverse as well:\n",
    "\n",
    "$$\\begin{bmatrix} 1 & 0 \\\\ 0 & 0 \\end{bmatrix}$$\n",
    "\n",
    "Why not? Well, just multiply it with a general $2 \\times 2$ matrix and see what happens to the value on the lower right corner of the result. Regardless of the values you put in the other matrix, you will always get a $0$ there, and it should be a $1$! Of course maybe this does not seem surprising, as this matrix has _\"lots of zeroes\"_. But this matrix has no zeros and it is also **not** invertible:\n",
    "\n",
    "$$\\begin{bmatrix} 1 & \\frac{2}{3} \\\\ 3 & 2 \\end{bmatrix}$$\n",
    "\n",
    "And this matrix has a couple of zeroes and is invertible:\n",
    "\n",
    "$$\\begin{bmatrix} 3 & 0 \\\\ 0 & \\pi \\end{bmatrix}$$\n",
    "\n",
    "We will leave the implementation of the inverse of a matrix for the section [Inverses revisited](#Inverses-revisited), after we have read about [determinants](#Determinant) and about [Gaussian elimination](#Gaussian elimination)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinant\n",
    "\n",
    "To determine if a matrix is invertible, we can check its _determinant_. Just like a number $x$ is invertible when $x \\neq 0$, we have that a matrix is invertible when its determinant is different from $0$, which we write $\\det{A} \\neq 0$ if we are talking about a matrix $A$ being invertible. But computing the determinant of a matrix is not a trivial task like looking at a number and deciding if it is $0$ or not. At least generally.\n",
    "\n",
    "In this section I will often include remarks about the [big O complexity](https://en.wikipedia.org/wiki/Big_O_notation) of some operations, in terms of the size $n$ of the matrix. If you are not familiar with this, do **not** worry as those will be just side comments to support my point that computing the determinant can be quite _hard_ in various senses.\n",
    "\n",
    "If your matrix is a diagonal matrix, then its determinant is just the product of all the entries in the diagonal. If we recall the notation $diag(a_1, a_2, \\cdots, a_n)$ for a diagonal matrix with elements $a_1, a_2, \\cdots, a_n$ in the diagonal, then\n",
    "\n",
    "$$\\det(diag(a_1, a_2, \\cdots, a_n)) = \\prod_{i=1}^n a_i$$\n",
    "\n",
    " > If this is the case, then computing the determinant is going to be $\\mathcal{O}(n)$ but you probably have to check if the matrix is diagonal before, which is done in $\\mathcal{O}(n^2)$, so everything together becomes $\\mathcal{O}(n^2)$ for a diagonal matrix.\n",
    " \n",
    "There are also two more types of matrices for which the determinant is easy to compute. These are the _upper triangular_ and _lower triangular_ matrices, which are \"matrices with lots of zeros\". But those zeroes cannot be scattered accross the matrix! A _lower triangular matrix_ $L$ is a matrix where all values above the diagonal are $0$, so it is a matrix of this form\n",
    "\n",
    "$$L = \\begin{bmatrix} a_{11} & 0 & 0 & \\cdots & 0 \\\\\n",
    "                        a_{21} & a_{22} & 0 &  &  \\\\\n",
    "                        \\vdots &     & a_{33} & \\ddots & \\vdots \\\\\n",
    "                              &       &      & \\ddots & 0 \\\\\n",
    "                        a_{n1} &  & \\cdots & & a_{nn}\\end{bmatrix}$$\n",
    "                        \n",
    "where the idea is that the interesting values are all in the _lower_ part of the matrix. Similarly, an _upper triangular matrix_ $U$ is a matrix that is $0$ everywhere below the diagonal. So, if this made sense to you then we can agree that a diagonal matrix is a matrix that is both upper triangular and lower triangular. For both these types of matrices the determinant is also computed by multiplying all the values in the diagonal. It is not obvious why this is the case right now but it will become obvious later on.\n",
    "\n",
    " > Finding out if a matrix is upper/lower triangular is $\\mathcal{O}(n^2)$ and then multiplying the values in the diagonal is $\\mathcal{O}(n)$ so we have again that computing the determinant is $\\mathcal{O}(n^2)$.\n",
    " \n",
    "The actual definition of determinant is, from the linear algebra standpoint, quite interesting. The determinant can even be motivated as a generalization of the length of a line, area of a rectangle and volume of a parallelepiped we are used to. In the end, the formula that comes out of this motivated definition is [#terrible](https://en.wikipedia.org/wiki/Determinant#n_%C3%97_n_matrices): both to understand and from the computational standpoint. So we will naturally skip all the subtleties there are to it and present a formula that works for a general matrix $A$.\n",
    "\n",
    "In order to do that, we need two final remarks:\n",
    " - The determinant of a $1 \\times 1$ matrix (that is, a scalar) is the number itself.\n",
    " - We need to introduce the minors $M_i(A)$ of a matrix $A$.\n",
    " \n",
    "If $A$ is an $n \\times n$ matrix, then its minor $M_{i,j}(A)$ is the $(n-1)\\times(n-1)$ matrix that you obtain by deleting the $i$-th row and $j$-th column of $A$. For example,\n",
    "\n",
    "$$A = \\begin{bmatrix} 1 & 2 \\\\ 3 & 4 \\end{bmatrix} \\implies M_{1,1}(A) = 4$$\n",
    "\n",
    "or for a bigger example,\n",
    "\n",
    "$$A = \\begin{bmatrix} 1 & 2 & 3 & 4 \\\\ 5 & 6 & 7 & 8 \\\\ 9 & 10 & 11 & 12 \\\\ 13 & 14 & 15 & 16 \\end{bmatrix} \\implies\n",
    "M_{3,2}(A) = \\begin{bmatrix} 1 & 3 & 4 \\\\ 5 & 7 & 8 \\\\ 13 & 15 & 16 \\end{bmatrix}$$\n",
    "\n",
    "With the minors defined, we can give a _recurrence_ formula (known as [Laplace's formula](https://en.wikipedia.org/wiki/Laplace_expansion)) for the determinant of $A$.\n",
    "\n",
    "$$\\det A = \\sum_{i=1}^n (-1)^{i+j}a_{i,j}\\det M_{i,j}(A) = \\sum_{j=1}^n (-1)^{i+j}a_{i,j}\\det M_{i,j}(A)$$\n",
    "\n",
    "Notice that for the first equality, $1 \\leq j \\leq n$ is fixed and for the second equality, $1 \\leq i \\leq n$ is fixed. In practice you can \"go through\" a row or \"go through\" a column and compute some minors along that row/column. If you pick a row/column with many zeroes, then most of the terms will be zero from the start.\n",
    "\n",
    " > It may not obvious at a first glance but this recurrence formula is $\\mathcal{O}(n!)$, which is one of the worst performances you could ask for in an algorithm... This formula has this performance because computing the determinant of a $n \\times n$ matrix involves computing $n$ determinants of $(n-1)\\times(n-1)$ matrices.\n",
    "\n",
    "This can be seen as the reason why computing the determinant of a lower/upper triangular matrix is just multiplying all the values in the diagonal. Try it out! Take a random lower triangular matrix and try using [Laplace's formula](https://en.wikipedia.org/wiki/Laplace_expansion) over the first row to compute the determinant. All terms other than the first are going to be zero. And then you have to compute the determinant of a minor, which will also be a lower triangular matrix! If you keep going, you end up multiplying all the values in the diagonal.\n",
    " \n",
    "Well, this was a long section! Now off to some coding! We got quite a lot done!\n",
    " - On lines 160-162 we defined a function `is_diagonal` that tells if a matrix is a diagonal matrix, and this function uses\n",
    "  - The function `is_lower_triangular` on lines 164-170\n",
    "  - The function `is_upper_triangular` on lines 172-178\n",
    " - On lines 180-187 we defined the function `minor(i,j)` that returns a minor of the matrix\n",
    " - On lines 189-204 we defined the determinant of a matrix given by Laplace's formula with the `determinant` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[12, 4]\n",
      " [0, 123]]\n",
      "False\n",
      "True\n",
      "[[1, 0, 0]\n",
      " [2, 3, 0]\n",
      " [4, 5, 6]]\n",
      "True\n",
      "False\n",
      "<Matrix (6x6)>\n",
      "True\n",
      "True\n",
      "True\n",
      "720\n",
      "[[1, 2, 3, 4]\n",
      " [5, 6, 7, 8]\n",
      " [9, 10, 11, 12]\n",
      " [13, 14, 15, 16]]\n",
      "[[1, 3, 4]\n",
      " [5, 7, 8]\n",
      " [13, 15, 16]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies a matrix with another matrix with a scalar\"\"\"\n",
    "        if not isinstance(other, (Matrix, int, float, complex)):\n",
    "            raise TypeError(\"Cannot multiply Matrix with {}\".format(type(other)))\n",
    "        elif isinstance(other, Matrix) and self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        \n",
    "        if isinstance(other, Matrix):\n",
    "            # this gives the shape of the result:\n",
    "            nrows = self.nrows()\n",
    "            ncols = other.ncols()\n",
    "            values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "            for r, row in enumerate(self.get_rows()):\n",
    "                for c, col in enumerate(other.get_cols()):\n",
    "                    values[r][c] = Matrix.dot_product(row, col)\n",
    "            return Matrix.from_values(values)\n",
    "        else:\n",
    "            return Matrix.from_values([\n",
    "                [value*other for value in row] for row in self.get_rows()\n",
    "            ])\n",
    "        \n",
    "    def __rmul__(self, other):\n",
    "        \"\"\"This is called when self was on the right-hand side of the * operator\n",
    "            as in 5*Matrix or 4.2*Matrix\"\"\"\n",
    "        return self.__mul__(other)\n",
    "    \n",
    "    def is_diagonal(self):\n",
    "        \"\"\"Returns true if the matrix is a diagonal matrix\"\"\"\n",
    "        return self.is_upper_triangular() and self.is_lower_triangular()\n",
    "    \n",
    "    def is_lower_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is a lower triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i+1, self.ncols()):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def is_upper_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is an upper triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def minor(self, ii, jj):\n",
    "        \"\"\"Returns the (ii,jj)th minor of the matrix\n",
    "        Indexes the matrix by 1 and not by 0\"\"\"\n",
    "        new_values = [\n",
    "            [self._values[i][j] for j in range(self.ncols()) if j+1 != jj]\n",
    "                                for i in range(self.nrows()) if i+1 != ii\n",
    "        ]\n",
    "        return Matrix.from_values(new_values)\n",
    "    \n",
    "    def determinant(self):\n",
    "        \"\"\"Returns the determinant of the matrix\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise TypeError(\"Cannot compute the determinant of a non-square matrix!\")\n",
    "        elif self.nrows() == 1:\n",
    "            return self._values[0][0]\n",
    "        else:\n",
    "            acc = 0\n",
    "            # use the formula with the minors over the first column\n",
    "            for i in range(1, self.ncols()+1):\n",
    "                # if this factor gives 0, skip it\n",
    "                if self._values[i-1][0] == 0:\n",
    "                    continue\n",
    "                minor = self.minor(i, 1)\n",
    "                acc += (-1)**(i+1)*self._values[i-1][0]*minor.determinant()\n",
    "            return acc\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\"\n",
    "    \n",
    "m = Matrix.from_values([[12,4],[0, 123]])\n",
    "print(m.full_str())\n",
    "print(m.is_lower_triangular())\n",
    "print(m.is_upper_triangular())\n",
    "\n",
    "m = Matrix.from_values([[1, 0, 0],[2, 3, 0], [4, 5, 6]])\n",
    "print(m.full_str())\n",
    "print(m.is_lower_triangular())\n",
    "print(m.is_upper_triangular())\n",
    "\n",
    "m = Matrix.from_diagonal([1,2,3,4,5,6])\n",
    "print(m)\n",
    "print(m.is_lower_triangular())\n",
    "print(m.is_upper_triangular())\n",
    "print(m.is_diagonal())\n",
    "print(m.determinant())\n",
    "\n",
    "m = Matrix.from_values([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])\n",
    "print(m.full_str())\n",
    "print(m.minor(3,2).full_str())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transpose\n",
    "\n",
    "There is another operation that is defined on matrices and it may seem weird. Not because it is complicated, but because it is not obvious why we would care. Either way, said operation is the _transpose_. The _transpose_ of a matrix $A$ that is $n \\times m$ is the matrix $A^T$ of shape $m \\times n$ satisfying the following equality:\n",
    "\n",
    "$$(A^T)_{i,j} = A_{j,i}$$\n",
    "\n",
    "which means that the element in the $i$-th row, $j$-th column of the transpose of $A$ is the element that was in the $j$-th row, $i$-th column of $A$. You can see this transpose operation as flipping the matrix along its diagonal. So, if the matrix $A$ is square, then $diag(A) = diag(A^T)$. As an example,\n",
    "\n",
    "$$A = \\begin{bmatrix} 1 & 2 & 3 & 4 \\\\ 5 & 6 & 7 & 8 \\\\ 9 & 10 & 11 & 12 \\\\ 13 & 14 & 15 & 16 \\end{bmatrix} \\implies A^T = \\begin{bmatrix} 1 & 5 & 9 & 13 \\\\ 2 & 6 & 10 & 14 \\\\ 3 & 7 & 11 & 15 \\\\ 4 & 8 & 12 & 16\\end{bmatrix}$$\n",
    "\n",
    "But $A$ does not need to be square; for example, the tranpose of a row vector is a column vector, and vice-versa:\n",
    "\n",
    "$$v = \\begin{bmatrix} 1 & 2 & 4 \\end{bmatrix} \\implies v^T \\begin{bmatrix} 1 \\\\ 2 \\\\ 4 \\end{bmatrix} \\implies (v^T)^T = \\begin{bmatrix} 1 & 2 & 4 \\end{bmatrix}$$\n",
    "\n",
    "From the previous example, and if you give it a little thought, it should become obvious that the transpose of the transpose always equals the original matrix. So $(A^T)^T = A$. But when is the transpose of a matrix equal to the matrix itself? When does $A^T = A$? Of course it can only happen for square matrices, otherwise $A^T$ and $A$ don't even have the same shape; for the cases where $A$ is square and its transpose is equal to itself, we say that $A$ is symmetric.\n",
    "\n",
    "$$A \\text{ is symmetric} \\implies A^T = A$$\n",
    "\n",
    "If instead we have that $A^T = -A$, then we say that $A$ is antisymmetric. One final remark: if $A$ is antisymmetric, then $A$ has a diagonal filled with zeroes.\n",
    "\n",
    "$$A \\text{ is antisymmetric} \\implies A^T = -A \\text{ and } diag(A) = [0, 0, \\cdots, 0]$$\n",
    "\n",
    "Now we will take care of implementing the transpose operation and an equality check for matrices.\n",
    " - In lines 95-107 we implemented the `__eq__` function that allows to check if two matrices are the same\n",
    " - In lines 174-176 we implemented the `__neg__` function, which is the one that allows one to write `-M` for a `Matrix` `M`\n",
    " - In lines 195-200 we implemented the `transpose` function that returns the transpose of the matrix\n",
    " - In lines 182-186 and 188-193 we implemented respectively the `is_symmetric` and `is_antisymmetric` functions that check if a matrix has the two properties mentioned\n",
    " \n",
    "(To see the line numbering, choose a code cell and then press \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n",
      "[[1, 2, 3]\n",
      " [4, 5, 6]]\n",
      "[[1, 4]\n",
      " [2, 5]\n",
      " [3, 6]]\n",
      "--------------------\n",
      "[[1, 2, 3]\n",
      " [2, 4, 5]\n",
      " [3, 5, 6]]\n",
      "True\n",
      "False\n",
      "--------------------\n",
      "[[0, 1, 2]\n",
      " [-1, 0, 3]\n",
      " [-2, -3, 0]]\n",
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __eq__(self, other):\n",
    "        \"\"\"Checks if two matrices are the same\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            return False\n",
    "        elif self.ncols() != self.ncols() or self.nrows() != other.nrows():\n",
    "            return False\n",
    "        these_cols = self.get_cols()\n",
    "        other_cols = other.get_cols()\n",
    "        for c1, c2 in zip(these_cols, other_cols):\n",
    "            for elem1, elem2 in zip(c1, c2):\n",
    "                if elem1 != elem2:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies a matrix with another matrix with a scalar\"\"\"\n",
    "        if not isinstance(other, (Matrix, int, float, complex)):\n",
    "            raise TypeError(\"Cannot multiply Matrix with {}\".format(type(other)))\n",
    "        elif isinstance(other, Matrix) and self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        \n",
    "        if isinstance(other, Matrix):\n",
    "            # this gives the shape of the result:\n",
    "            nrows = self.nrows()\n",
    "            ncols = other.ncols()\n",
    "            values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "            for r, row in enumerate(self.get_rows()):\n",
    "                for c, col in enumerate(other.get_cols()):\n",
    "                    values[r][c] = Matrix.dot_product(row, col)\n",
    "            return Matrix.from_values(values)\n",
    "        else:\n",
    "            return Matrix.from_values([\n",
    "                [value*other for value in row] for row in self.get_rows()\n",
    "            ])\n",
    "        \n",
    "    def __rmul__(self, other):\n",
    "        \"\"\"This is called when self was on the right-hand side of the * operator\n",
    "            as in 5*Matrix or 4.2*Matrix\"\"\"\n",
    "        return self.__mul__(other)\n",
    "    \n",
    "    def __neg__(self):\n",
    "        \"\"\"For a matrix M returns -M\"\"\"\n",
    "        return -1*self\n",
    "    \n",
    "    def is_diagonal(self):\n",
    "        \"\"\"Returns true if the matrix is a diagonal matrix\"\"\"\n",
    "        return self.is_upper_triangular() and self.is_lower_triangular()\n",
    "    \n",
    "    def is_symmetric(self):\n",
    "        \"\"\"Returns true if the matrix is a symmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        return self == self.transpose()\n",
    "    \n",
    "    def is_antisymmetric(self):\n",
    "        \"\"\"Returns true if the matrix is an antisymmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        # need to define __neg__ for this\n",
    "        return self == -self.transpose()\n",
    "    \n",
    "    def transpose(self):\n",
    "        \"\"\"Returns the transpose of this matrix\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in col] for col in self.get_cols()\n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def is_lower_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is a lower triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i+1, self.ncols()):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def is_upper_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is an upper triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def minor(self, ii, jj):\n",
    "        \"\"\"Returns the (ii,jj)th minor of the matrix\n",
    "        Indexes the matrix by 1 and not by 0\"\"\"\n",
    "        new_values = [\n",
    "            [self._values[i][j] for j in range(self.ncols()) if j+1 != jj]\n",
    "                                for i in range(self.nrows()) if i+1 != ii\n",
    "        ]\n",
    "        return Matrix.from_values(new_values)\n",
    "    \n",
    "    def determinant(self):\n",
    "        \"\"\"Returns the determinant of the matrix\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise TypeError(\"Cannot compute the determinant of a non-square matrix!\")\n",
    "        elif self.nrows() == 1:\n",
    "            return self._values[0][0]\n",
    "        else:\n",
    "            acc = 0\n",
    "            # use the formula with the minors over the first column\n",
    "            for i in range(1, self.ncols()+1):\n",
    "                # if this factor gives 0, skip it\n",
    "                if self._values[i-1][0] == 0:\n",
    "                    continue\n",
    "                minor = self.minor(i, 1)\n",
    "                acc += (-1)**(i+1)*self._values[i-1][0]*minor.determinant()\n",
    "            return acc\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\"\n",
    "    \n",
    "v1 = [[1, 2, 3], [4, 5, 6]]\n",
    "v2 = [[1, 2, 3], [0, 5, 6]]\n",
    "m1 = Matrix.from_values(v1)\n",
    "m2 = Matrix.from_values(v2)\n",
    "m3 = Matrix.from_values(v2)\n",
    "print(m1 == m2)\n",
    "print(m2 == m3)\n",
    "print(m1.full_str())\n",
    "print(m1.transpose().full_str())\n",
    "print(\"-\"*20)\n",
    "v4 = [[1, 2, 3], [2, 4, 5], [3, 5, 6]]\n",
    "m4 = Matrix.from_values(v4)\n",
    "print(m4.full_str())\n",
    "print(m4.is_symmetric())\n",
    "print(m4.is_antisymmetric())\n",
    "print(\"-\"*20)\n",
    "v5 = [[0, 1, 2], [-1, 0, 3], [-2, -3, 0]]\n",
    "m5 = Matrix.from_values(v5)\n",
    "print(m5.full_str())\n",
    "print(m5.is_symmetric())\n",
    "print(m5.is_antisymmetric())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear systems\n",
    "\n",
    "One of the many reasons why matrices are so useful is because they make it easier to see if a linear system has a solution; and when it does, matrices also make it easier for us to solve it.\n",
    "\n",
    "Take this example linear system of equations:\n",
    "$$\\begin{cases}\n",
    "3x + 2y = 0 \\\\\n",
    "-3y + z = 3 \\\\\n",
    "4x - z = 12\\end{cases}$$\n",
    "\n",
    "Now we will encode this linear system as $Ax = b$ where $x$ and $b$ are column vectors and $A$ is a $3 \\times 3$ matrix:\n",
    "$$\\begin{cases}\n",
    "3x + 2y = 0\\\\ \n",
    "-3y + z = 3\\\\\n",
    "4x - z = 12\\end{cases} \\iff \\begin{bmatrix} 3 & 2 & 0 \\\\ 0 & -3 & 1 \\\\ 4 & 0 & -1 \\end{bmatrix}\\begin{bmatrix}x \\\\ y \\\\ z\\end{bmatrix} = \\begin{bmatrix}0 \\\\ 3 \\\\ 12 \\end{bmatrix} \\iff Ax = b $$\n",
    "Take a moment to realize that the matrix multiplication I wrote, with $Ax$, actually gives back the linear system. Now imagine that $A$ was actually an invertible matrix! Then, $Ax = b \\iff A^{-1}Ax = b \\iff x = A^{-1}b$ and thus you only need to do the matrix multiplication $A^{-1}b$ to find the values of your variables.\n",
    "\n",
    "Not only that, but we know that $A^{-1}$ only exists if $\\det A \\neq 0$, so now we actually have a way of checking if a linear system of $n$ equations and $n$ unknowns will be solvable or not.\n",
    "\n",
    "### Gaussian elimination\n",
    "\n",
    "Right now I just wrote\n",
    " > you _only_ need to do the matrix multiplication $A^{-1}b$ to find the values of your variables\n",
    "\n",
    "It just so happens that you **also** have to find the inverse of $A$! Which can be quite difficult... So in practice we should use other methods! The method of [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) is surely the one that is taught the most.\n",
    "\n",
    "The method, at first sight, might not be trivial, but it grows on you as you use it more and more. The idea is quite simple. If you have two equalities, namely\n",
    "\n",
    "$$\\begin{cases}\n",
    "x = y \\\\\n",
    "a = b\\end{cases}$$\n",
    "\n",
    "then it is also true that\n",
    "\n",
    "$$ a + w x = b + w y $$\n",
    "\n",
    "Please take your time to digest that new equality I wrote. If you think straight, it is not that surprising. Now the objective of Gaussian elimination is to consecutively write \"new equalities\" by taking the equalities we already have and summing multiples of those to the other equalities. If we make our decisions in a smart way, we can simplify our equations in the process.\n",
    "\n",
    "For the sake of simplicity we will assume that our gaussian elimination method will only be asked to solve systems that have one and only one solution. If a user wants to know if this is the case, the user can always check the determinant of the matrix that encodes the linear system first. If the determinant is **not** zero, then the system has one and only one solution. We will also assume the matrix $Ax = b$ in $A$ is square. This is not a strict requirement for the system $Ax = b$ to have one and only one solution, but it makes our lives much easier.\n",
    "\n",
    "Other than that, let us just say that if `M` is a `Matrix` object then `M.linear_solve(b)` attempts to solve the linear system `Mx = b`. Now let us write some code!\n",
    " - In lines 42-48 we implemented the `copy` function that creates a [deep copy](https://stackoverflow.com/questions/184710/what-is-the-difference-between-a-deep-copy-and-a-shallow-copy) of a Matrix.\n",
    " - In lines 252-309 we implement our longest function so far, `linear_solve(b)` that tries to solve the linear system `self*x = b` with respect to `x`.\n",
    " \n",
    "(To see the line numbering, choose a code cell and then press \"L\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "False\n",
      "##########\n",
      "[[10.0]\n",
      " [3.0]]\n",
      "[[10]\n",
      " [3]]\n",
      "##########\n",
      "[[10.0]\n",
      " [3.0]\n",
      " [-3.0000000000000004]\n",
      " [0.3333333333333428]]\n",
      "[[10]\n",
      " [3]\n",
      " [-3]\n",
      " [0.3333333333333333]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def copy(self):\n",
    "        \"\"\"Returns a copy M of the matrix self, such that\n",
    "        `M == self` is true but `M is self` is False\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in row] for row in self.get_rows() \n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __eq__(self, other):\n",
    "        \"\"\"Checks if two matrices are the same\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            return False\n",
    "        elif self.ncols() != self.ncols() or self.nrows() != other.nrows():\n",
    "            return False\n",
    "        these_cols = self.get_cols()\n",
    "        other_cols = other.get_cols()\n",
    "        for c1, c2 in zip(these_cols, other_cols):\n",
    "            for elem1, elem2 in zip(c1, c2):\n",
    "                if elem1 != elem2:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies a matrix with another matrix with a scalar\"\"\"\n",
    "        if not isinstance(other, (Matrix, int, float, complex)):\n",
    "            raise TypeError(\"Cannot multiply Matrix with {}\".format(type(other)))\n",
    "        elif isinstance(other, Matrix) and self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        \n",
    "        if isinstance(other, Matrix):\n",
    "            # this gives the shape of the result:\n",
    "            nrows = self.nrows()\n",
    "            ncols = other.ncols()\n",
    "            values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "            for r, row in enumerate(self.get_rows()):\n",
    "                for c, col in enumerate(other.get_cols()):\n",
    "                    values[r][c] = Matrix.dot_product(row, col)\n",
    "            return Matrix.from_values(values)\n",
    "        else:\n",
    "            return Matrix.from_values([\n",
    "                [value*other for value in row] for row in self.get_rows()\n",
    "            ])\n",
    "        \n",
    "    def __rmul__(self, other):\n",
    "        \"\"\"This is called when self was on the right-hand side of the * operator\n",
    "            as in 5*Matrix or 4.2*Matrix\"\"\"\n",
    "        return self.__mul__(other)\n",
    "    \n",
    "    def __neg__(self):\n",
    "        \"\"\"For a matrix M returns -M\"\"\"\n",
    "        return -1*self\n",
    "    \n",
    "    def is_diagonal(self):\n",
    "        \"\"\"Returns true if the matrix is a diagonal matrix\"\"\"\n",
    "        return self.is_upper_triangular() and self.is_lower_triangular()\n",
    "    \n",
    "    def is_symmetric(self):\n",
    "        \"\"\"Returns true if the matrix is a symmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        return self == self.transpose()\n",
    "    \n",
    "    def is_antisymmetric(self):\n",
    "        \"\"\"Returns true if the matrix is an antisymmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        # need to define __neg__ for this\n",
    "        return self == -self.transpose()\n",
    "    \n",
    "    def transpose(self):\n",
    "        \"\"\"Returns the transpose of this matrix\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in col] for col in self.get_cols()\n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def is_lower_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is a lower triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i+1, self.ncols()):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def is_upper_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is an upper triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def minor(self, ii, jj):\n",
    "        \"\"\"Returns the (ii,jj)th minor of the matrix\n",
    "        Indexes the matrix by 1 and not by 0\"\"\"\n",
    "        new_values = [\n",
    "            [self._values[i][j] for j in range(self.ncols()) if j+1 != jj]\n",
    "                                for i in range(self.nrows()) if i+1 != ii\n",
    "        ]\n",
    "        return Matrix.from_values(new_values)\n",
    "    \n",
    "    def determinant(self):\n",
    "        \"\"\"Returns the determinant of the matrix\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise TypeError(\"Cannot compute the determinant of a non-square matrix!\")\n",
    "        elif self.nrows() == 1:\n",
    "            return self._values[0][0]\n",
    "        else:\n",
    "            acc = 0\n",
    "            # use the formula with the minors over the first column\n",
    "            for i in range(1, self.ncols()+1):\n",
    "                # if this factor gives 0, skip it\n",
    "                if self._values[i-1][0] == 0:\n",
    "                    continue\n",
    "                minor = self.minor(i, 1)\n",
    "                acc += (-1)**(i+1)*self._values[i-1][0]*minor.determinant()\n",
    "            return acc\n",
    "        \n",
    "    def linear_solve(self, b):\n",
    "        \"\"\"Attempts to solve the square linear system `self*X = b`\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise ValueError(\"Uh oh, linear_solve expected a square matrix!\")\n",
    "        n = self.nrows()\n",
    "        # we will perform gaussian elimination in a copy of self\n",
    "        # extended with a copy of the target b\n",
    "        rows_self = self.get_rows()\n",
    "        rows_b = b.get_rows()\n",
    "        if len(rows_self) != len(rows_b):\n",
    "            raise ValueError(\"Uh oh, linear_solve expected the target b \" + \\\n",
    "                                \"to have {} rows\".format(len(rows_self)))\n",
    "        extended_rows = [r1 + r2 for r1, r2 in zip(rows_self, rows_b)]\n",
    "        M = Matrix.from_values(extended_rows)\n",
    "        # first turn the first part of M into an upper triangular matrix\n",
    "        for i in range(n):\n",
    "            # try to place a non-zero element at the (i,i) position\n",
    "            idx = i\n",
    "            column = M.get_col(i+1)\n",
    "            while idx < n and column[idx] == 0:\n",
    "                idx += 1\n",
    "            else:\n",
    "                if idx >= n:\n",
    "                    raise ValueError(\"Uh oh, linear_solve expected a solvable system\")\n",
    "            # do we have to swap rows `i` and `idx` to get a non-zero diagonal element?\n",
    "            if idx != i:\n",
    "                rowi = M.get_row(i+1)[:]\n",
    "                rowidx = M.get_row(idx+1)[:]\n",
    "                values = M._values\n",
    "                values[i] = rowidx\n",
    "                values[idx] = rowi\n",
    "                M = Matrix.from_values(values)\n",
    "            # go down the column, zero-ing everything below element (i,i)\n",
    "            rowi = M.get_row(i+1)\n",
    "            values = M._values\n",
    "            for j in range(i+1, n):\n",
    "                mult = -values[j][i]/values[i][i]\n",
    "                values[j] = [mult*rowi_elem + val \\\n",
    "                             for val, rowi_elem in zip(values[j], rowi)]\n",
    "            M = Matrix.from_values(values)\n",
    "        # then turn the first part of M into the id matrix\n",
    "        values = M._values\n",
    "        for i in range(n-1, -1, -1):\n",
    "            # first rescale the i-th row to get a 1 in the diagonal\n",
    "            diag = values[i][i]\n",
    "            values[i] = [elem/diag for elem in values[i]]\n",
    "            # go up the column, zero-ing everything else above (i,i)\n",
    "            rowi = values[i][:]\n",
    "            for j in range(i):\n",
    "                mult = -values[j][i]\n",
    "                values[j] = [mult*rowi_elem + val \\\n",
    "                            for val, rowi_elem in zip(values[j], rowi)]\n",
    "        # at this point `values` holds the result\n",
    "        # extract the second part\n",
    "        result_values = [\n",
    "            row[n:] for row in values\n",
    "        ]\n",
    "        return Matrix.from_values(result_values)\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\"\n",
    "\n",
    "v1 = [[1, 2, 3], [4, 5, 6]]\n",
    "m1 = Matrix.from_values(v1)\n",
    "m = m1.copy()\n",
    "print(m == m1) # should be True\n",
    "print(m is m1)\n",
    "###\n",
    "print(\"#\"*10)\n",
    "###\n",
    "vA = [[1, 2], [3, 4]]\n",
    "A = Matrix.from_values(vA)\n",
    "vb = [[10], [3]]\n",
    "b = Matrix.from_values(vb)\n",
    "x = A.linear_solve(b)\n",
    "print((A*x).full_str())\n",
    "print(b.full_str())\n",
    "###\n",
    "print(\"#\"*10)\n",
    "###\n",
    "vA = [[0, 2, 0, 10], [0, 4, 0, 25], [1, 0, -3, 0.2], [15, 15, 15, 0]]\n",
    "A = Matrix.from_values(vA)\n",
    "vb = [[10], [3], [-3], [1/3]]\n",
    "b = Matrix.from_values(vb)\n",
    "x = A.linear_solve(b)\n",
    "print((A*x).full_str())\n",
    "print(b.full_str())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inverses revisited\n",
    "\n",
    "It is time to write a function that returns the inverse of a matrix! Your first task is to go check the function `linear_solve` and try to realize what the nature of the argument `b` can be.\n",
    "\n",
    "We wrote the line `b.get_rows()` so clearly `b` must be a `Matrix`, and also the length of `b.get_rows()` should match the length of `self.get_rows()`, meaning that the number of rows of `b` should match the number of rows of the matrix. But other than that, we didn't specify how many columns `b` can/must have! This can be quite handy! In fact, if we want to solve the systems $Ax_1 = b_1$ and $Ax_2 = b_2$ where the matrix $A$ is the same, we can effectively build a matrix $B = [b_1\\ b_2]$ with two columns, $b_1$ and $b_2$ and then call `A.linear_solve(B)`. The result will then be a matrix $X$ with two columns, $X = [x_1\\ x_2]$. The code just works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b1\n",
      "[[-17.0]\n",
      " [5.0]]\n",
      "b2\n",
      "[[-14.0]\n",
      " [4.0]]\n",
      "B\n",
      "[[-17.0, -14.0]\n",
      " [5.0, 4.0]]\n"
     ]
    }
   ],
   "source": [
    "vA = [[0, 2], [1, 4]]\n",
    "A = Matrix.from_values(vA)\n",
    "# b1\n",
    "vb1 = [[10], [3]]\n",
    "b1 = Matrix.from_values(vb1)\n",
    "x1 = A.linear_solve(b1)\n",
    "print(\"b1\")\n",
    "print(x1.full_str())\n",
    "# b2\n",
    "vb2 = [[8], [2]]\n",
    "b2 = Matrix.from_values(vb2)\n",
    "x2 = A.linear_solve(b2)\n",
    "print(\"b2\")\n",
    "print(x2.full_str())\n",
    "# B\n",
    "vB = [[10, 8], [3, 2]]\n",
    "B = Matrix.from_values(vB)\n",
    "X = A.linear_solve(B)\n",
    "print(\"B\")\n",
    "print(X.full_str())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is just great! We can use this extra functionality we got for free to build the inverse function. In fact, think of it this way, if $A$ is a $n\\times n$ matrix with an inverse $B$, then $B$ has $n$ columns that we will call $b_1, b_2, \\cdots, b_n$. Also, let us agree to call $e_i$ to the column vector that is $0$ everywhere _except_ in position $i$, where it is a $1$ (these are usually called the [vectors of the standard basis](https://en.wikipedia.org/wiki/Standard_basis)).\n",
    "\n",
    "So, for example, if $n = 3$ then\n",
    "\n",
    "$$e_1 = \\begin{bmatrix} 1 \\\\ 0 \\\\ 0 \\end{bmatrix} \\quad e_2 = \\begin{bmatrix} 0 \\\\ 1 \\\\ 0 \\end{bmatrix} \\quad e_3 = \\begin{bmatrix} 0 \\\\ 0 \\\\ 1 \\end{bmatrix}$$\n",
    "\n",
    "I introduced these vectors $e_i$ because now, if we really understand the way matrix multiplication works and what it means for $B$ to be the inverse of $A$, we should have $Ab_i = e_i$. Do you agree? Maybe write down the two multiplications with\n",
    "\n",
    "$$A = \\begin{bmatrix} 1 & 2 \\\\ 0 & 1 \\end{bmatrix}, \\quad B = \\begin{bmatrix}1 & -2 \\\\ 0 & 1\\end{bmatrix}$$\n",
    "\n",
    "So, basically we just want to solve the $n$ linear systems $Ab_i = e_i$. But because of the way our function `linear_solve` works, we can just solve one big system. For that we need to put all of the $e_i$ together, but that creates the identity matrix $I_n$! So we can finally define our `inverse` function, that we did on lines 311-323."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.855, 0.74, 0.25, 0.05]\n",
      " [2.5, -1.0, 0.0, 0.0]\n",
      " [-0.645, 0.26, -0.25, 0.016666666666666666]\n",
      " [-0.4, 0.2, -0.0, -0.0]]\n",
      "[[1.0, 0.0, 0.0, 0.0]\n",
      " [0.0, 1.0, 0.0, 0.0]\n",
      " [5.551115123125783e-17, -2.7755575615628914e-17, 1.0, 0.0]\n",
      " [0.0, 0.0, 0.0, 1.0]]\n"
     ]
    }
   ],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def copy(self):\n",
    "        \"\"\"Returns a copy M of the matrix self, such that\n",
    "        `M == self` is true but `M is self` is False\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in row] for row in self.get_rows() \n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __eq__(self, other):\n",
    "        \"\"\"Checks if two matrices are the same\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            return False\n",
    "        elif self.ncols() != self.ncols() or self.nrows() != other.nrows():\n",
    "            return False\n",
    "        these_cols = self.get_cols()\n",
    "        other_cols = other.get_cols()\n",
    "        for c1, c2 in zip(these_cols, other_cols):\n",
    "            for elem1, elem2 in zip(c1, c2):\n",
    "                if elem1 != elem2:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies a matrix with another matrix with a scalar\"\"\"\n",
    "        if not isinstance(other, (Matrix, int, float, complex)):\n",
    "            raise TypeError(\"Cannot multiply Matrix with {}\".format(type(other)))\n",
    "        elif isinstance(other, Matrix) and self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        \n",
    "        if isinstance(other, Matrix):\n",
    "            # this gives the shape of the result:\n",
    "            nrows = self.nrows()\n",
    "            ncols = other.ncols()\n",
    "            values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "            for r, row in enumerate(self.get_rows()):\n",
    "                for c, col in enumerate(other.get_cols()):\n",
    "                    values[r][c] = Matrix.dot_product(row, col)\n",
    "            return Matrix.from_values(values)\n",
    "        else:\n",
    "            return Matrix.from_values([\n",
    "                [value*other for value in row] for row in self.get_rows()\n",
    "            ])\n",
    "        \n",
    "    def __rmul__(self, other):\n",
    "        \"\"\"This is called when self was on the right-hand side of the * operator\n",
    "            as in 5*Matrix or 4.2*Matrix\"\"\"\n",
    "        return self.__mul__(other)\n",
    "    \n",
    "    def __neg__(self):\n",
    "        \"\"\"For a matrix M returns -M\"\"\"\n",
    "        return -1*self\n",
    "    \n",
    "    def is_diagonal(self):\n",
    "        \"\"\"Returns true if the matrix is a diagonal matrix\"\"\"\n",
    "        return self.is_upper_triangular() and self.is_lower_triangular()\n",
    "    \n",
    "    def is_symmetric(self):\n",
    "        \"\"\"Returns true if the matrix is a symmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        return self == self.transpose()\n",
    "    \n",
    "    def is_antisymmetric(self):\n",
    "        \"\"\"Returns true if the matrix is an antisymmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        # need to define __neg__ for this\n",
    "        return self == -self.transpose()\n",
    "    \n",
    "    def transpose(self):\n",
    "        \"\"\"Returns the transpose of this matrix\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in col] for col in self.get_cols()\n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def is_lower_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is a lower triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i+1, self.ncols()):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def is_upper_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is an upper triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def minor(self, ii, jj):\n",
    "        \"\"\"Returns the (ii,jj)th minor of the matrix\n",
    "        Indexes the matrix by 1 and not by 0\"\"\"\n",
    "        new_values = [\n",
    "            [self._values[i][j] for j in range(self.ncols()) if j+1 != jj]\n",
    "                                for i in range(self.nrows()) if i+1 != ii\n",
    "        ]\n",
    "        return Matrix.from_values(new_values)\n",
    "    \n",
    "    def determinant(self):\n",
    "        \"\"\"Returns the determinant of the matrix\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise TypeError(\"Cannot compute the determinant of a non-square matrix!\")\n",
    "        elif self.nrows() == 1:\n",
    "            return self._values[0][0]\n",
    "        else:\n",
    "            acc = 0\n",
    "            # use the formula with the minors over the first column\n",
    "            for i in range(1, self.ncols()+1):\n",
    "                # if this factor gives 0, skip it\n",
    "                if self._values[i-1][0] == 0:\n",
    "                    continue\n",
    "                minor = self.minor(i, 1)\n",
    "                acc += (-1)**(i+1)*self._values[i-1][0]*minor.determinant()\n",
    "            return acc\n",
    "        \n",
    "    def linear_solve(self, b):\n",
    "        \"\"\"Attempts to solve the square linear system `self*X = b`\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise ValueError(\"Uh oh, linear_solve expected a square matrix!\")\n",
    "        n = self.nrows()\n",
    "        # we will perform gaussian elimination in a copy of self\n",
    "        # extended with a copy of the target b\n",
    "        rows_self = self.get_rows()\n",
    "        rows_b = b.get_rows()\n",
    "        if len(rows_self) != len(rows_b):\n",
    "            raise ValueError(\"Uh oh, linear_solve expected the target b \" + \\\n",
    "                                \"to have {} rows\".format(len(rows_self)))\n",
    "        extended_rows = [r1 + r2 for r1, r2 in zip(rows_self, rows_b)]\n",
    "        M = Matrix.from_values(extended_rows)\n",
    "        # first turn the first part of M into an upper triangular matrix\n",
    "        for i in range(n):\n",
    "            # try to place a non-zero element at the (i,i) position\n",
    "            idx = i\n",
    "            column = M.get_col(i+1)\n",
    "            while idx < n and column[idx] == 0:\n",
    "                idx += 1\n",
    "            else:\n",
    "                if idx >= n:\n",
    "                    raise ValueError(\"Uh oh, linear_solve expected a solvable system\")\n",
    "            # do we have to swap rows `i` and `idx` to get a non-zero diagonal element?\n",
    "            if idx != i:\n",
    "                rowi = M.get_row(i+1)[:]\n",
    "                rowidx = M.get_row(idx+1)[:]\n",
    "                values = M._values\n",
    "                values[i] = rowidx\n",
    "                values[idx] = rowi\n",
    "                M = Matrix.from_values(values)\n",
    "            # go down the column, zero-ing everything below element (i,i)\n",
    "            rowi = M.get_row(i+1)\n",
    "            values = M._values\n",
    "            for j in range(i+1, n):\n",
    "                mult = -values[j][i]/values[i][i]\n",
    "                values[j] = [mult*rowi_elem + val \\\n",
    "                             for val, rowi_elem in zip(values[j], rowi)]\n",
    "            M = Matrix.from_values(values)\n",
    "        # then turn the first part of M into the id matrix\n",
    "        values = M._values\n",
    "        for i in range(n-1, -1, -1):\n",
    "            # first rescale the i-th row to get a 1 in the diagonal\n",
    "            diag = values[i][i]\n",
    "            values[i] = [elem/diag for elem in values[i]]\n",
    "            # go up the column, zero-ing everything else above (i,i)\n",
    "            rowi = values[i][:]\n",
    "            for j in range(i):\n",
    "                mult = -values[j][i]\n",
    "                values[j] = [mult*rowi_elem + val \\\n",
    "                            for val, rowi_elem in zip(values[j], rowi)]\n",
    "        # at this point `values` holds the result\n",
    "        # extract the second part\n",
    "        result_values = [\n",
    "            row[n:] for row in values\n",
    "        ]\n",
    "        return Matrix.from_values(result_values)\n",
    "    \n",
    "    def inverse(self):\n",
    "        \"\"\"Tries to return the inverse of the matrix.\n",
    "        If no inverse exists, an error is raised\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise ValueError(\"Cannot invert non-square matrix\")\n",
    "        n = self.nrows()\n",
    "        Id = Matrix.eye(n)\n",
    "        try:\n",
    "            inv = self.linear_solve(Id)\n",
    "        # catch the error to change the message\n",
    "        except ValueError:\n",
    "            raise ValueError(\"Cannot invert matrix\")\n",
    "        return inv\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\"\n",
    "\n",
    "vA = [[0, 2, 0, 10], [0, 4, 0, 25], [1, 0, -3, 0.2], [15, 15, 15, 0]]\n",
    "A = Matrix.from_values(vA)\n",
    "inv = A.inverse()\n",
    "print(inv.full_str())\n",
    "# notice the annoying floating point errors\n",
    "print((A*inv).full_str())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Matrix implementation\n",
    "\n",
    "For your convenience, below you can find the final version of the Matrix object you built. Notice that the [numpy](https://www.numpy.org/) module is a module that does everything we just implemented, and more! And in a much more efficient way! But according to the **No free lunch theorem**, if we used that from the get-go we'd be losing out on something... Which is the deeper understanding you get if you truly try to implement these functions on your own and understand how they work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Matrix(object):\n",
    "    def __init__(self, rows, columns, values=None):\n",
    "        \"\"\"Initialize a matrix with its number of rows and columns\"\"\"\n",
    "        # the _ is a standard that says these variables\n",
    "        # are not to be accessed directly from outside\n",
    "        self._rows = rows\n",
    "        self._cols = columns\n",
    "        \n",
    "        if values is None:\n",
    "            self._values = [[0 for j in range(columns)] for i in range(rows)]\n",
    "        else:\n",
    "            # We were provided with some initial values\n",
    "            #  copy them one by one, as doing self._values = values\n",
    "            #  has some unintended side-effects\n",
    "            # Ensure the sizes match\n",
    "            if len(values) != rows:\n",
    "                raise ValueError(\"Size mismatch rows {} != {}\".format(rows, len(values)))\n",
    "            self._values = []\n",
    "            for row in values:\n",
    "                matrix_row = []\n",
    "                if len(row) != columns:\n",
    "                    raise ValueError(\"Size mismatch columns {} != {}\".format(columns, len(row)))\n",
    "                for value in row:\n",
    "                    matrix_row.append(value)\n",
    "                self._values.append(matrix_row)\n",
    "                \n",
    "    def from_values(values):\n",
    "        \"\"\"Takes the given list of values and converts to matrix;\n",
    "            Assumes the list of lists will be well-shaped\"\"\"\n",
    "        rows = len(values)\n",
    "        columns = len(values[0])\n",
    "        return Matrix(rows, columns, values)\n",
    "    \n",
    "    def from_diagonal(diagonal):\n",
    "        \"\"\"Takes a list with the values of the diagonal of a matrix;\n",
    "            Returns the corresponding diagonal matrix\"\"\"\n",
    "        values = [[0 for i in range(len(diagonal))] for j in range(len(diagonal))]\n",
    "        for i, value in enumerate(diagonal):\n",
    "            values[i][i] = value\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def copy(self):\n",
    "        \"\"\"Returns a copy M of the matrix self, such that\n",
    "        `M == self` is true but `M is self` is False\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in row] for row in self.get_rows() \n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def eye(n):\n",
    "        \"\"\"Returns the identity matrix of shape NxN\"\"\"\n",
    "        diag = [1 for i in range(n)]\n",
    "        return Matrix.from_diagonal(diag)\n",
    "    \n",
    "    def shape(self):\n",
    "        \"\"\"Returns the shape of the matrix as (nrows, ncols)\"\"\"\n",
    "        return (self.nrows(), self.ncols())\n",
    "        \n",
    "    def nrows(self):\n",
    "        \"\"\"Return the number of rows of the matrix\"\"\"\n",
    "        return self._rows\n",
    "    \n",
    "    def get_row(self, index):\n",
    "        \"\"\"Return the specified row\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Row index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Row index {} cannot be < 1\".format(index))\n",
    "        elif index > self.nrows():\n",
    "            raise IndexError(\"Row index {} cannot be > {}\".format(index, self.nrows()))\n",
    "        # make a deep copy of the row\n",
    "        row = []\n",
    "        for value in self._values[index-1]:\n",
    "            row.append(value)\n",
    "        return row\n",
    "    \n",
    "    def get_rows(self):\n",
    "        \"\"\"Return a list of all rows\"\"\"\n",
    "        return [self.get_row(r) for r in range(1, self.nrows()+1)]\n",
    "    \n",
    "    def ncols(self):\n",
    "        \"\"\"Return the number of columns of the matrix\"\"\"\n",
    "        return self._cols\n",
    "    \n",
    "    def get_col(self, index):\n",
    "        \"\"\"Return the specified column\"\"\"\n",
    "        if not isinstance(index, int):\n",
    "            raise TypeError(\"Column index should be an integer\")\n",
    "        elif index < 1:\n",
    "            raise IndexError(\"Column index {} cannot be < 1\".format(index))\n",
    "        elif index > self.ncols():\n",
    "            raise IndexError(\"Column index {} cannot be > {}\".format(index, self.ncols()))\n",
    "        # make a deep copy of the column\n",
    "        column = []\n",
    "        for i in range(self.nrows()):\n",
    "            column.append(self._values[i][index-1])\n",
    "        return column\n",
    "    \n",
    "    def get_cols(self):\n",
    "        \"\"\"Return a list of all columns\"\"\"\n",
    "        return [self.get_col(c) for c in range(1, self.ncols()+1)]\n",
    "    \n",
    "    def __eq__(self, other):\n",
    "        \"\"\"Checks if two matrices are the same\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            return False\n",
    "        elif self.ncols() != self.ncols() or self.nrows() != other.nrows():\n",
    "            return False\n",
    "        these_cols = self.get_cols()\n",
    "        other_cols = other.get_cols()\n",
    "        for c1, c2 in zip(these_cols, other_cols):\n",
    "            for elem1, elem2 in zip(c1, c2):\n",
    "                if elem1 != elem2:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        \"\"\"Adds two matrices together\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            # the \"other\" is not a matrix, what should I do..?\n",
    "            raise TypeError(\"Can only add a Matrix with another Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        # we have two matrices and they have the same shape, great!\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA + valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        \"\"\"Subtracts two matrices; notice the operation is self - other\"\"\"\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise TypeError(\"Can only subtract a Matrix from a Matrix\")\n",
    "        elif self.shape() != other.shape():\n",
    "            raise ValueError(\"Matrices should have the same shape: {} != {}\".format(self.shape(), other.shape()))\n",
    "        newMatrix = []\n",
    "        for rowA, rowB in zip(self._values, other._values):\n",
    "            newRow = []\n",
    "            for valA, valB in zip(rowA, rowB):\n",
    "                newRow.append(valA - valB)\n",
    "            newMatrix.append(newRow)\n",
    "        return Matrix.from_values(newMatrix)\n",
    "    \n",
    "    def dot_product(u, v):\n",
    "        \"\"\"Computes the dot product between vectors u and v\n",
    "        Given as lists\"\"\"\n",
    "        acc = 0\n",
    "        for ui, vi in zip(u, v):\n",
    "            acc += ui*vi\n",
    "        return acc\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"Multiplies a matrix with another matrix with a scalar\"\"\"\n",
    "        if not isinstance(other, (Matrix, int, float, complex)):\n",
    "            raise TypeError(\"Cannot multiply Matrix with {}\".format(type(other)))\n",
    "        elif isinstance(other, Matrix) and self.ncols() != other.nrows():\n",
    "            raise Error(\"Number of columns on left {} should match \" + \\\n",
    "                            \"number of rows on right {}\".format(self.ncols(), other.nrows()))\n",
    "        \n",
    "        if isinstance(other, Matrix):\n",
    "            # this gives the shape of the result:\n",
    "            nrows = self.nrows()\n",
    "            ncols = other.ncols()\n",
    "            values = [[0 for c in range(ncols)] for r in range(nrows)]\n",
    "            for r, row in enumerate(self.get_rows()):\n",
    "                for c, col in enumerate(other.get_cols()):\n",
    "                    values[r][c] = Matrix.dot_product(row, col)\n",
    "            return Matrix.from_values(values)\n",
    "        else:\n",
    "            return Matrix.from_values([\n",
    "                [value*other for value in row] for row in self.get_rows()\n",
    "            ])\n",
    "        \n",
    "    def __rmul__(self, other):\n",
    "        \"\"\"This is called when self was on the right-hand side of the * operator\n",
    "            as in 5*Matrix or 4.2*Matrix\"\"\"\n",
    "        return self.__mul__(other)\n",
    "    \n",
    "    def __neg__(self):\n",
    "        \"\"\"For a matrix M returns -M\"\"\"\n",
    "        return -1*self\n",
    "    \n",
    "    def is_diagonal(self):\n",
    "        \"\"\"Returns true if the matrix is a diagonal matrix\"\"\"\n",
    "        return self.is_upper_triangular() and self.is_lower_triangular()\n",
    "    \n",
    "    def is_symmetric(self):\n",
    "        \"\"\"Returns true if the matrix is a symmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        return self == self.transpose()\n",
    "    \n",
    "    def is_antisymmetric(self):\n",
    "        \"\"\"Returns true if the matrix is an antisymmetric matrix\"\"\"\n",
    "        if self.ncols() != self.nrows():\n",
    "            return False\n",
    "        # need to define __neg__ for this\n",
    "        return self == -self.transpose()\n",
    "    \n",
    "    def transpose(self):\n",
    "        \"\"\"Returns the transpose of this matrix\"\"\"\n",
    "        values = [\n",
    "            [elem for elem in col] for col in self.get_cols()\n",
    "        ]\n",
    "        return Matrix.from_values(values)\n",
    "    \n",
    "    def is_lower_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is a lower triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i+1, self.ncols()):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def is_upper_triangular(self):\n",
    "        \"\"\"Returns true if the matrix is an upper triangular matrix\"\"\"\n",
    "        for i in range(self.nrows()):\n",
    "            for j in range(i):\n",
    "                if self._values[i][j] != 0:\n",
    "                    return False\n",
    "        return True\n",
    "    \n",
    "    def minor(self, ii, jj):\n",
    "        \"\"\"Returns the (ii,jj)th minor of the matrix\n",
    "        Indexes the matrix by 1 and not by 0\"\"\"\n",
    "        new_values = [\n",
    "            [self._values[i][j] for j in range(self.ncols()) if j+1 != jj]\n",
    "                                for i in range(self.nrows()) if i+1 != ii\n",
    "        ]\n",
    "        return Matrix.from_values(new_values)\n",
    "    \n",
    "    def determinant(self):\n",
    "        \"\"\"Returns the determinant of the matrix\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise TypeError(\"Cannot compute the determinant of a non-square matrix!\")\n",
    "        elif self.nrows() == 1:\n",
    "            return self._values[0][0]\n",
    "        else:\n",
    "            acc = 0\n",
    "            # use the formula with the minors over the first column\n",
    "            for i in range(1, self.ncols()+1):\n",
    "                # if this factor gives 0, skip it\n",
    "                if self._values[i-1][0] == 0:\n",
    "                    continue\n",
    "                minor = self.minor(i, 1)\n",
    "                acc += (-1)**(i+1)*self._values[i-1][0]*minor.determinant()\n",
    "            return acc\n",
    "        \n",
    "    def linear_solve(self, b):\n",
    "        \"\"\"Attempts to solve the square linear system `self*X = b`\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise ValueError(\"Uh oh, linear_solve expected a square matrix!\")\n",
    "        n = self.nrows()\n",
    "        # we will perform gaussian elimination in a copy of self\n",
    "        # extended with a copy of the target b\n",
    "        rows_self = self.get_rows()\n",
    "        rows_b = b.get_rows()\n",
    "        if len(rows_self) != len(rows_b):\n",
    "            raise ValueError(\"Uh oh, linear_solve expected the target b \" + \\\n",
    "                                \"to have {} rows\".format(len(rows_self)))\n",
    "        extended_rows = [r1 + r2 for r1, r2 in zip(rows_self, rows_b)]\n",
    "        M = Matrix.from_values(extended_rows)\n",
    "        # first turn the first part of M into an upper triangular matrix\n",
    "        for i in range(n):\n",
    "            # try to place a non-zero element at the (i,i) position\n",
    "            idx = i\n",
    "            column = M.get_col(i+1)\n",
    "            while idx < n and column[idx] == 0:\n",
    "                idx += 1\n",
    "            else:\n",
    "                if idx >= n:\n",
    "                    raise ValueError(\"Uh oh, linear_solve expected a solvable system\")\n",
    "            # do we have to swap rows `i` and `idx` to get a non-zero diagonal element?\n",
    "            if idx != i:\n",
    "                rowi = M.get_row(i+1)[:]\n",
    "                rowidx = M.get_row(idx+1)[:]\n",
    "                values = M._values\n",
    "                values[i] = rowidx\n",
    "                values[idx] = rowi\n",
    "                M = Matrix.from_values(values)\n",
    "            # go down the column, zero-ing everything below element (i,i)\n",
    "            rowi = M.get_row(i+1)\n",
    "            values = M._values\n",
    "            for j in range(i+1, n):\n",
    "                mult = -values[j][i]/values[i][i]\n",
    "                values[j] = [mult*rowi_elem + val \\\n",
    "                             for val, rowi_elem in zip(values[j], rowi)]\n",
    "            M = Matrix.from_values(values)\n",
    "        # then turn the first part of M into the id matrix\n",
    "        values = M._values\n",
    "        for i in range(n-1, -1, -1):\n",
    "            # first rescale the i-th row to get a 1 in the diagonal\n",
    "            diag = values[i][i]\n",
    "            values[i] = [elem/diag for elem in values[i]]\n",
    "            # go up the column, zero-ing everything else above (i,i)\n",
    "            rowi = values[i][:]\n",
    "            for j in range(i):\n",
    "                mult = -values[j][i]\n",
    "                values[j] = [mult*rowi_elem + val \\\n",
    "                            for val, rowi_elem in zip(values[j], rowi)]\n",
    "        # at this point `values` holds the result\n",
    "        # extract the second part\n",
    "        result_values = [\n",
    "            row[n:] for row in values\n",
    "        ]\n",
    "        return Matrix.from_values(result_values)\n",
    "    \n",
    "    def inverse(self):\n",
    "        \"\"\"Tries to return the inverse of the matrix.\n",
    "        If no inverse exists, an error is raised\"\"\"\n",
    "        if self.nrows() != self.ncols():\n",
    "            raise ValueError(\"Cannot invert non-square matrix\")\n",
    "        n = self.nrows()\n",
    "        Id = Matrix.eye(n)\n",
    "        try:\n",
    "            inv = self.linear_solve(Id)\n",
    "        # catch the error to change the message\n",
    "        except ValueError:\n",
    "            raise ValueError(\"Cannot invert matrix\")\n",
    "        return inv\n",
    "    \n",
    "    def __str__(self):\n",
    "        \"\"\"String representation of the matrix for printing\"\"\"\n",
    "        return \"<Matrix ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def __repr__(self):\n",
    "        \"\"\"Representation of the matrix for printing inside containers\"\"\"\n",
    "        return \"<M ({}x{})>\".format(self.nrows(), self.ncols())\n",
    "    \n",
    "    def full_str(self):\n",
    "        \"\"\"Returns a full string representation of the values of the matrix\"\"\"\n",
    "        return \"[\"+\"\\n \".join(list(map(str, self.get_rows())))+\"]\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Contacts\n",
    "\n",
    "If you enjoyed this guide and/or it was useful, consider leaving a star in the [Virgilio repository](https://github.com/virgili0/Virgilio) and sharing it with your friends!\n",
    "\n",
    "This was brought to you by the editor of the [Mathspp Blog](https://mathspp.com/blog) (or [Mathspp on Facebook](https://www.facebook.com/mathspp/)), [RojerGS](https://github.com/RojerGS). Write to me at mathsppblog @ gmail.com."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
